#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <unistd.h>
#include <math.h>
#include <string.h>
#include <pthread.h>
#include <sys/time.h>

#define DEBUG 0
#define MAXTHREADS 32

typedef struct
{
	unsigned int id;
	
	double *matrix;
	double *vector;
	double *result;
	int dim;
	int size;

} thread_data;

void usage(void)
{
	printf("Syntax Error!\n");
	printf("Usage: ./mvmult <matrix dimension (50..N)> <iterations (1..N)>\n");
	exit(-1);
}

enum times {START, END, DIFF};
double calcRuntime(struct timeval *t)
{
	/* calc time difference */
	t[DIFF].tv_sec = t[END].tv_sec - t[START].tv_sec;
	if (t[START].tv_usec > t[END].tv_usec)
	{
		t[DIFF].tv_usec = (1e6 - t[END].tv_usec) + t[START].tv_usec;
	}
	else
	{
		t[DIFF].tv_usec = t[END].tv_usec - t[START].tv_usec;
	}
	
	return ((double)t[DIFF].tv_sec + (double)t[DIFF].tv_usec/1e6);
}

void mvmult(double *m, double *v, double *r, int dim, int worksize)
{
	unsigned int i, j;

	//printf("dim=%i worksize=%i\n", dim, worksize);
	for (i=0; i<worksize; i++)
	{
		for (j=0; j<dim; j++)
		{
			*(r + i) += *(m + (i*dim) + j) * *(v + j);
			//printf("(%.1lf * %.1lf) + ", *(m + (i*dim) + j), *(v + j) );
		}
		//printf("\n");
	}
	//printf("thread end\n\n\n");
	/* Use usleep() when calling the threads to get a useful dbg-output here */
}

void* worker(void* data)
{
	thread_data *td;
	td = (thread_data*)data;
	//printf("#%i: this is thread %i\n", td->id, td->id);
	//printf("#%i: size=%i\n", td->id, td->size);
	//printf("#%i:  dim=%i\n", td->id, td->dim);
	mvmult(td->matrix, td->vector, td->result, td->dim, td->size);

	return NULL;
}

int main(int argc, char *argv[])
{
	if (argc != 3) usage();
	
	unsigned int i, j;
	struct timeval t[3];
	double runtime;
	double seq_runtime = 1;

	int nThreads;
	int nDim = atoi(argv[1]);
	int iterations = atoi(argv[2]);
	if (nDim <= 50 || iterations <= 0) usage();
	
	/* allocate space */
	double *matrix, *vector, *result;
	matrix = (double*)malloc(nDim*nDim*sizeof(double));
	vector = (double*)malloc(nDim*sizeof(double));
	result = (double*)malloc(nDim*sizeof(double));

	thread_data td[MAXTHREADS];	
	pthread_t threads[MAXTHREADS];

	for (nThreads=1; nThreads<=MAXTHREADS; nThreads<<=1)
	{
		/* initialize */
		for (i=0; i<nDim; i++)
		{
			*(vector + i) = (rand() % 1000) * 0.1f;
			//*(vector + i) = i;
			*(result + i) = 0;
			for (j=0; j<nDim; j++)
			{
				*(matrix + (i*nDim) + j) = (rand() % 1000) * 0.1f;
				//*(matrix + (i*nDim) + j) = j;
			}
		}

#if DEBUG == 1
		printf("Vector: \n");
		for (i=0; i<nDim; i++)
		{
			printf("%3.1lf\n", *(vector + i));
		}
	
		printf("Matrix: \n");
		for (i=0; i<nDim; i++)
		{
			for (j=0; j<nDim; j++)
			{
				printf("%3.1lf ", *(matrix + (i*nDim) + j));
			}
			printf("\n");
		}
	
		printf("Result: \n");
		for (i=0; i<nDim; i++)
		{
			printf("%3.1lf\n", *(result + i));
		}
#endif

		/* distribute work */
		int workPerThread, workRest;
		workPerThread = nDim / nThreads;
		workRest = nDim % nThreads;
		//printf("# threads=%i wpt=%i, wr=%i\n", nThreads, workPerThread, workRest);

		for (i=0; i<nThreads; i++)
		{
			td[i].id = i;
	
			td[i].dim = nDim;
			td[i].size = workPerThread;
			td[i].vector = vector;
			td[i].matrix = matrix + (i*workPerThread*nDim);
			td[i].result = result + (i*workPerThread);
		}
	
		runtime = 0.0;
		for (int it=0; it<iterations; it++)
		{
			//double percent = (double)it/(double)iterations *100+.1;
			//printf("\r# progress: %.1lf %% done", percent);
			memset(result, 0, nDim*sizeof(double));
			gettimeofday(&t[0], NULL);
#if 1
			/* start workers*/
			for (i=0; i<nThreads; i++)
			{
				pthread_create(&threads[i], NULL, worker, (void*)(&td[i]));
				//usleep(1000);
			}
		
			/* compute the rest */
			if (workRest != 0)
			{
				int m_offset = nThreads*nDim*workPerThread;
				int r_offset = nThreads*workPerThread;

				//printf("# rest present\n");
				mvmult(matrix+m_offset, vector, result+r_offset, nDim, workRest);
			}

			/* join workers */
			for (i=0; i<nThreads; i++)
			{
				pthread_join(threads[i], NULL);
			}
#else
			/* single-thread call for reference:*/
			mvmult(matrix, vector, result, nDim, nDim);
#endif
			gettimeofday(&t[1], NULL);

#if DEBUG == 1
			printf("\n# Result: \n");
			for (i=0; i<nDim; i++)
			{
				printf("# %4.1lf\n", *(result + i));
			}
			printf("\n");
#endif

			runtime += calcRuntime(t);
		}
		//printf("\n");

		runtime /= iterations;
		
		if (nThreads == 1)
		{
			seq_runtime = runtime;
			printf("# seq=%lf\n", seq_runtime);
		}
		
		double speedup = seq_runtime/runtime;
		
		printf("%i %i %7.6lf %3.2lf %i\n", nThreads, nDim, runtime, speedup, iterations);	
	}
	printf("# threads, dim, runtime, speedup, iterations\n\n");
	
	free(matrix);
	free(vector);
	free(result);
	
	return EXIT_SUCCESS;
}

